Age and Behavior-Dependent Differential miRNAs Expression in the Hypopharyngeal Glands of Honeybees (Apis mellifera L.)

Simple Summary The hypopharyngeal glands (HPGs) are a pair of aciniform glands that are located in the frontal area of the heads of worker bees (Apis mellifera L.) that exhibit age and behavior-dependent development. Little is known about whether/how miRNAs regulate the HPGs development. In this study, small RNA sequencing was employed to analyze the miRNA profiles of HPGs in newly-emerged bees (NEB), nurse bees (NB), and forager bees (FB). We found that there were a total of 31 known miRNAs differentially expressed among the three stages, which might have regulatory roles in the growth and development, protein synthesis, and carbohydrate and energy metabolism in the HPGs. Additionally, the downregulation of ame-miR-184-3p and ame-miR-252a-5p in nurse bees may be involved in royal jelly secretion, while the lower expression of ame-miR-11-3p and ame-miR-281-3p in forager bees are responsible for honey processing. Abstract This study aims to investigate the expression differences of miRNAs in the hypopharyngeal glands (HPGs) of honeybees at three developmental stages and to explore their regulation functions in the HPGs development. Small RNA sequencing was employed to analyze the miRNA profiles of HPGs in newly-emerged bees (NEB), nurse bees (NB), and forager bees (FB). Results showed that a total of 153 known miRNAs were found in the three stages, and ame-miR-276-3p, ame-miR-375-3p, ame-miR-14-3p, ame-miR-275-3p, and ame-miR-3477-5p were the top five most abundant ones. Furthermore, the expression of 11 miRNAs, 17 miRNAs, and 18 miRNAs were significantly different in NB vs. FB comparison, NB vs. NEB comparison, and in FB vs. NEB comparison, respectively, of which ame-miR-184-3p and ame-miR-252a-5p were downregulated in NB compared with that in both the FB and NEB, while ame-miR-11-3p, ame-miR-281-3p, and ame-miR-31a-5p had lower expression levels in FB compared with that in both the NB and NEB. Bioinformatic analysis showed that the potential target genes of the differentially expressed miRNAs (DEMs) were mainly enriched in several key signaling pathways, including mTOR signaling pathway, MAPK signaling pathway-fly, FoxO signaling pathway, Hippo signaling pathway-fly. Overall, our study characterized the miRNA profiles in the HPGs of honeybees at three different developmental stages and provided a basis for further study of the roles of miRNAs in HPGs development.


Introduction
Honeybees (Apis mellifera L.) are eusocial insects, and a typical colony consists of three castes, including queen, drones, and worker bees. Amongst those, worker bees are the predominant ones that exhibit a complex age-dependent division of labor [1]. In general, within the first two weeks after emergence, the worker bees mainly perform tasks inside the hive, such as brood and queen rearing (nurse bees). After that, most of the

Small RNA Sequencing and Analysis
A total of three micrograms of RNA per sample was used to prepare the sequencing sample. Sequencing libraries were generated using NEBNext R Ultra R small RNA Sample Library Prep Kit for Illumina R (New England Biolabs, Beijing, China) following manufacturer's recommendations, and index codes were added to attribute sequences to each sample. In total, nine sequencing libraries were constructed, and each group (NEB, NB, and FB) contained three libraries. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v4-cBot-HS (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 4000 platform and 1 × 50 bp single-end reads were generated (Illumina). After quality control, the filter ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), other ncRNA. and repeat sequences were removed based on high-quality clean reads mapped to Silva database, GtRNAdb database, Rfam database, and Repbase database using Bowtie software [31]. Then, the remaining reads were used to detect known miRNA by comparing with known miRNAs in A. mellifera from miRbase, while RNAfold software [32] was used for novel miRNA secondary structure prediction.

Differential Expression Analysis
The miRNA expression levels were calculated by transcript per million (TPM). Differential expression analysis between the following comparisons: (1) NB vs. FB, (2) NB vs. NEB, and (3) FB vs. NEB, were calculated by the DESeq R package [33]. The miRNAs with a p-value < 0.01 and |log2-fold change| > 1 were assigned as significantly differentially expressed. The 3 -untranslated region (3 -UTR) sequences from the A. mellifera genome were used to predict the target genes of the differentially expressed miRNAs (DEMs) with miRanda [34] and TargetScan [35] software. The predicted target genes with miranda_energy < −10 or TargetScan_score ≥ 50 were obtained. We extracted the overlapping target genes identified using the two pieces of software. Gene ontology (GO) enrichment analysis of the predicted target genes was implemented by the topGO R package [36], while we used KOBAS software [37] to test the statistical enrichment of predicted target genes in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Small RNA Sequencing Analysis
Using the Illumina Hiseq 4000 platform, nine sRNA libraries were subjected to deep sRNA sequencing. After removing the useless sequences, 3.27 to 6.85 million clean reads were obtained (Table S1). The length distribution of the small RNAs for all three groups peaked at 22 nt, with 16.87% to 31.36% of total clean reads ( Figure 1).

Small RNA Sequencing Analysis
Using the Illumina Hiseq 4000 platform, nine sRNA libraries were subjected to deep sRNA sequencing. After removing the useless sequences, 3.27 to 6.85 million clean reads were obtained (Table S1). The length distribution of the small RNAs for all three groups peaked at 22 nt, with 16.87% to 31.36% of total clean reads ( Figure 1).
After mapping to Silva database, GtRNAdb database, Rfam database, and Repbase database, the clean reads were divided into different categories, and 60.71-72.30% belonged to unique rRNAs, 5.85-15.54% were unique tRNAs, 0.87-2.83% were unique snRNAs, and 0.50-1.47% were unique snoRNAs. On the other hand, the remaining 17.56-20.97% unannotated reads containing miRNAs were used to identify the known miRNAs and novel miRNAs.

Identification of Known and Novel miRNAs
We compared the unannotated sequences with known mature miRNAs of A. mellifera in the miRBase database to identify the known miRNAs. In total, 140, 130, and 129 known miRNAs were found in NB HPGs, FB HPGs, and NEB HPGs, respectively. Among them, 114 were shared between three groups ( Figure 2A, Table S2). Ten miRNAs with the highest expression were selected and listed in Table 2, and five of them (ame-miR-276-3p, ame-miR-375-3p, ame-miR-14-3p, ame-miR-275-3p, and ame-miR-3477-5p) were After mapping to Silva database, GtRNAdb database, Rfam database, and Repbase database, the clean reads were divided into different categories, and 60.71-72.30% belonged to unique rRNAs, 5.85-15.54% were unique tRNAs, 0.87-2.83% were unique snRNAs, and 0.50-1.47% were unique snoRNAs. On the other hand, the remaining 17.56-20.97% unannotated reads containing miRNAs were used to identify the known miRNAs and novel miRNAs.
In total, 26, 18, and 21 novel miRNAs were predicted in NB HPGs, FB HPGs, and NEB HPGs, respectively. Among them, six novel miRNAs were shared in the three groups, while seven, nine, and ten novel miRNAs were common in NB vs. FB, FB vs. NEB, and NB vs. NEB comparisons, respectively. Otherwise, fifteen, eight, and eight novel miRNAs were only expressed in NB HPGs, FB HPGs, and NEB HPGs, respectively ( Figure  2B, Table S3).

Differential Expression Analysis of miRNAs
As shown in Table 3, there were 31 miRNAs differentially expressed among the NEB, NB, and FB libraries (p < 0.01 and |log2-fold change| > 1). Among these, 11 miRNAs exhibited differential expression in the NB vs. FB comparison (five upregulated and six downregulated), 17 miRNAs in NB vs. NEB comparison (five upregulated and twelve downregulated), and 18 miRNAs in FB vs. NEB comparison (six upregulated and twelve downregulated). In addition, one miRNA, ame-miR-31a-5p, was common in the three comparisons; five miRNAs were common between NB vs. FB and NB vs. NEB comparisons, five miRNAs between NB vs. NEB and FB vs. NEB comparisons, and six miRNAs between FB vs. NEB and NB vs. NEB comparisons ( Figure 3, Table 3). Some differentially expressed miRNAs (DEMs) called our attention: ame-miR-3477-5p and ame-  In total, 26, 18, and 21 novel miRNAs were predicted in NB HPGs, FB HPGs, and NEB HPGs, respectively. Among them, six novel miRNAs were shared in the three groups, while seven, nine, and ten novel miRNAs were common in NB vs. FB, FB vs. NEB, and NB vs. NEB comparisons, respectively. Otherwise, fifteen, eight, and eight novel miRNAs were only expressed in NB HPGs, FB HPGs, and NEB HPGs, respectively ( Figure 2B, Table S3).

Prediction of the Target Genes of DEMs and Enrichment Analysis
Using miRanda and TargetScan softwares, we predicted the target genes of the 31 DEMs among the NEB, NB, and FB libraries. In total, 1180, 1531, and 1868 candidate target genes were obtained for 11 DEMs confirmed in the NB vs. FB comparison (Table S4), 17 DEMs confirmed in the NB vs. NEB comparison (Table S5), and 18 DEMs confirmed in the FB vs. NEB comparison (Table S6), respectively. The predicted target genes of these DEMs obtained in the three comparisons were enriched in thousands of GO terms. Nucleus, plasma membrane, protein binding, integral component of plasma membrane, and regulation of transcription, DNA-templated were the top five most enriched ones ( Figure 4, Tables S7-S9). The predicted target genes were mapped to 217 KEGG pathways, and the mTOR signaling pathway, MAPK signaling pathway-fly, FoxO signaling pathway, Hippo signaling pathway-fly, and Autophagy-animal were the top five most abundant ones ( Figure 5, Tables S10-S12).

Prediction of the Target Genes of DEMs and Enrichment Analysis
Using miRanda and TargetScan softwares, we predicted the target genes of the 31 DEMs among the NEB, NB, and FB libraries. In total, 1180, 1531, and 1868 candidate target genes were obtained for 11 DEMs confirmed in the NB vs. FB comparison (Table  S4), 17 DEMs confirmed in the NB vs. NEB comparison (Table S5), and 18 DEMs confirmed in the FB vs. NEB comparison (Table S6), respectively. The predicted target genes of these DEMs obtained in the three comparisons were enriched in thousands of GO terms. Nucleus, plasma membrane, protein binding, integral component of plasma membrane, and regulation of transcription, DNA-templated were the top five most enriched ones (Figure 4, Tables S7-S9). The predicted target genes were mapped to 217 KEGG pathways, and the mTOR signaling pathway, MAPK signaling pathway-fly, FoxO signaling pathway, Hippo signaling pathway-fly, and Autophagy-animal were the top five most abundant ones ( Figure 5, Tables S10-S12).

Discussion
In this study, sRNA-seq was employed to explore the miRNA profiles in the hypopharyngeal glands of honeybees at three developmental stages (newly-emerged bees, nurse bees, and forager bees). In total, 153 known miRNAs were obtained, and 114 of them were shared in the three stages. Differential expression analysis showed that 11, 17, and 18 miRNAs were significantly differentially expressed between NB vs. FB comparison, FB vs. NEB comparison, and NB vs. NEB comparison, respectively ( Figure 3, Table 3). Validation of ten DEMs using qPCR indicated that the results of qPCR and sRNA-seq showed a similar trend, suggesting the reliability of sRNA-seq data [39]. The potential target genes of DEMs were also predicted, which mainly participated in diverse signaling pathways, including mTOR signaling pathway, MAPK signaling pathway-fly,

Discussion
In this study, sRNA-seq was employed to explore the miRNA profiles in the hypopharyngeal glands of honeybees at three developmental stages (newly-emerged bees, nurse bees, and forager bees). In total, 153 known miRNAs were obtained, and 114 of them were shared in the three stages. Differential expression analysis showed that 11, 17, and 18 miRNAs were significantly differentially expressed between NB vs. FB comparison, FB vs. NEB comparison, and NB vs. NEB comparison, respectively ( Figure 3, Table 3). Validation of ten DEMs using qPCR indicated that the results of qPCR and sRNA-seq showed a similar trend, suggesting the reliability of sRNA-seq data [39]. The potential target genes of DEMs were also predicted, which mainly participated in diverse signaling pathways, including mTOR signaling pathway, MAPK signaling pathway-fly, FoxO signal-ing pathway, and Hippo signaling pathway-fly ( Figure 5, Tables S10-S12), which might play important roles in regulating HPGs development in honeybees.
Among the 153 known miRNAs obtained in the present study, ame-miR-276-3p, ame-miR-375-3p, ame-miR-14-3p, ame-miR-275-3p, and ame-miR-3477-5p were the top five most abundant miRNAs that commonly existed in NB, FB, and NEB libraries ( Table 2). A recent study demonstrated that ame-miR-14 was enriched in the queen ovary and played an important role in regulating egg-laying via modulating ecdysone receptor in honeybee queens [26]. The miRNA ame-mir-276 is highly expressed in the small-type Kenyon cells of the mushroom bodies and is involved with the development of related neural function in honeybees [40]. Therefore, the highly expressed miRNAs in our study might play indispensable roles in the HPGs development in honeybees. Furthermore, our study also found that the expression level of ame-miR-3477-5p exhibited significant differences in NB vs. FB comparison or in NB vs. NEB comparison, and ame-miR-276-3p had higher expression in NEB than that in NB. Previous studies reported that there were striking differences in the brain or head miRNAs expression between nurse and forager [21,22,41,42]. Hence, development stages have obvious impacts on the expression levels of miRNAs in honeybees [43,44].
There were a total of 31 differentially expressed miRNAs confirmed among the three different developmental stages, and several of them interested us enormously. For instance, ame-miR-184-3p and ame-miR-252a-5p had lower expression in NB than that in both the FB and NEB ( Table 3). The miRNA ame-miR-184-3p is found to be involved in royal jelly secretion by enhancing the expression of the target genes, insulin-like receptor and cyclin dependent kinase 12, and inhibition of ame-miR-184-3p expression also can activate the mTOR signaling pathway [45]. Accordingly, many potential target genes of ame-miR-184-3p obtained in the present study were enriched for this key signaling pathway (Table 4). In addition, many predicted target genes of ame-miR-184-3p and ame-miR-252a-5p were also responsible for protein synthesis and energy metabolism (Table 4), which is similar to a previous study declaring that the upregulated genes in the HPGs of NB were notably enriched in ribosome and aminoacyl-tRNA biosynthesis [15]. It is well known that the HPGs of nurse bees have strong activity to secrete RJ protein [13,46], and HPGs of high-RJproducing nurse bees exhibit stronger energy replenishment than those of Italian bees [47]. Hence, the lower expression of ame-miR-184-3p and ame-miR-252a-5p in nurse bees was closely associated with RJ secretion.
On the other hand, three miRNAs (ame-miR-11-3p, ame-miR-281-3p, and ame-miR-31a-5p) showed lower expression levels in FB compared with that in both the NB and NEB (Table 3). Similarly, a previous study found that ame-miR-31a had higher expression level in brains of NB than in brains of FB in both the typical colonies and single-cohort colonies, which was considered as an important regulator of the behavioral transition in worker bees [48]. The miRNAs can regulate various biological processes by suppressing their target gene expression [49]. Our present study found that several potential targets of ame-miR-11-3p (glucosidase 2 subunit beta and hexosaminidase D) and ame-miR-281-3p (myogenesis-regulating glycosidase) encoded the enzymes involved in carbohydrate digestion. Glucosidase can catalyze polysaccharides into glucose [50]. Hexosaminidase D participated in glycan degradation, and myogenesis-regulating glycosidase was responsible for starch and sucrose metabolism ( Table 4). The downregulation of ame-miR-11-3p and ame-miR-281-3p in HPGs of FB are likely involved with the conversion of nectar into honey [10,51]. Otherwise, some candidate target genes of ame-miR-11-3p, ame-miR-31a-5p, and ame-miR-281-3p were enriched in apoptosis-fly pathway (Table 4). Previous studies demonstrated that the HPGs begin to degenerate in worker bees more than 15 days old (forager bees) [52,53].
Bioinformatic analysis showed that mTOR signaling pathway, MAPK signaling pathwayfly, FoxO signaling pathway, Hippo signaling pathway-fly, and Autophagy-animal were the top five most enriched pathways in the potential target genes of the DEMs in the three comparisons: in NB vs. FB comparison, in NB vs. NEB comparison, and in FB vs. NEB comparison, which suggested that these pathways might play critical roles in hypopharyngeal glands development in honeybees. Among these, mTOR signaling has been implicated in the modulation of ribosome oogenesis and protein synthesis in cells [54,55]. Hippo signaling pathway is involved with organ size control and tissue regeneration via cell proliferation and apoptosis [56,57]. FoxO signaling pathway is responsible for carbohydrate and energy metabolism [58,59]. MAPK signaling pathway plays a significant role in regulating division of labor, caste differentiation, and queen development in honeybees [60][61][62]. These findings indicate that the DEMs might have regulatory roles in the growth and development, protein synthesis, and carbohydrate and energy metabolism in the HPGs of honeybees.

Conclusions
We first compared the miRNA profiles of hypopharyngeal glands in honeybees at different developmental stages (newly-emerged bees, nurse bees, and forager bees). Bioinformatic analysis showed that the differentially expressed miRNAs were involved in important biological processes related to growth and development, protein synthesis, and carbohydrate and energy metabolism in the hypopharyngeal glands. Additionally, we found that the downregulation of ame-miR-184-3p and ame-miR-252a-5p in nurse bees may be involved in royal jelly secretion, while the lower expression of ame-miR-11-3p and ame-miR-281-3p in forager bees are responsible for honey processing. These findings should provide a basis for further study of the roles of miRNAs in hypopharyngeal glands development.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/insects12090764/s1, Table S1: Details of the sRNA-seq data; Table S2: The known miRNAs detected in NB, FB, and NEB libraries; Table S3: The novel miRNAs detected in NB, FB, and NEB libraries; Table S4: Predicted target genes of differentially expressed miRNAs confirmed in the NB vs. FB comparison; Table S5: Predicted target genes of differentially expressed miRNAs confirmed in the NB vs. NEB comparison; Table S6: Predicted target genes of differentially expressed miRNAs confirmed in the FB vs. NEB comparison; Table S7: GO enrichment analysis of predicted target genes of differentially expressed miRNAs confirmed in the NB vs. FB comparison; Table S8: GO enrichment analysis of predicted target genes of differentially expressed miRNAs confirmed in the NB vs. NEB comparison; Table S9: GO enrichment analysis of predicted target genes of differentially expressed miRNAs confirmed in the FB vs. NEB comparison; Table S10: KEGG pathways analysis of predicted target genes of differentially expressed miRNAs confirmed in the NB vs. FB comparison; Table S11: KEGG pathways analysis of predicted target genes of differentially expressed miRNAs confirmed in the NB vs. NEB comparison; Table S12: KEGG pathways analysis of predicted target genes of differentially expressed miRNAs confirmed in the FB vs. NEB comparison.

Institutional Review Board Statement: Not applicable.
Data Availability Statement: Upon acceptance, the data used in this study will be available in the Supplementary Data File. The sequencing data are available in the SRA database (PRJNA754270) of the NCBI system.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, or the writing of the manuscript.